setwd(/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/20201207)
library(data.table)
sel=list.files(pattern="BrainSeqPhaseII")[-13]
filer=read.table("53K.add.GWAS.txt",head=T,sep="\t")

anno=fread("BrainSeqPhaseII_snp_annotation.txt",head=T,sep="\t")
anno=data.frame(snp=anno$snp,chr_hg19=anno$chr_hg19,pos_hg19=anno$pos_hg19)
anno$id=paste(anno$chr_hg19,anno$pos_hg19,sep=":")
anno=anno[anno$id %in% as.character(filer$id),]
for(i in 2:12){
file=fread(sel[i],head=T,sep="\t")
file=merge(file,anno,by="snp")
group.libd=paste(gsub(".txt","",unlist(strsplit(sel[i],"_"))[4]),gsub(".txt","",unlist(strsplit(sel[i],"_"))[5]),sep="_")
tmp=data.frame(snp=file$snp,id=file$id,group.LIBD=rep(group.libd,dim(file)[1]))
result=rbind(result,tmp)
}
result$id2=paste(result$id,result$group.LIBD,sep="_")
result = result[!duplicated(result$id2),]
result=result[,-4]

str1=unique(as.character(result$id))

rt=result[1,]
rt=rt[-1,]
for(j in 1:length(str1)){
tmp=result[result$id==str1[j],]
tmp1=tmp[1,]
tmp1$group.LIBD=paste(tmp$group.LIBD,collapse=";")
rt=rbind(rt,tmp1)
}
write.csv(rt,"53K.LIBD.eQTL.csv",quote=F,row.names=F)
